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ABSTRACT 

We present a data analysis pipeline for CMB polarization experiments, running from 
multi-frequency maps to the power spectra. We focus mainly on component separation 
and, for the first time, we work out the covariance matrix accounting for errors associ- 
ated to the separation itself. This allows us to propagate such errors and evaluate their 
contributions to the uncertainties on the final products. The pipeline is optimized for 
intermediate and small scales, but could be easily extended to lower multipoles. We 
exploit realistic simulations of the sky, tailored for the Planck mission. The compo- 
nent separation is achieved by exploiting the Correlated Component Analysis in the 
harmonic domain, t hat we demonstrate to be superior to the real-space application 
(|Bonaldi et al.ll2006l ). We present two techniques to estimate the uncertainties on the 
spectral parameters of the separated components. The component separation errors 
are then propagated by means of Monte Carlo simulations to obtain the corresponding 
contributions to uncertainties on the component maps and on the CMB power spec- 
tra. For the Planck polarization case they are found to be subdominant compared 
to noise. 



1 INTRODUCTION 



The Cosmic Microwave Background (CMB) radiation is a 
gold-mine of cosmological information. While cosmology is 
entering its precision era, the target of CMB experiments is 
shifting towards weak signals. The tiny amount of polariza- 
tion associated with the CMB anisotropy is undoubtedly one 
of the most intriguing - and challenging - measurements of 
this kind. Several experimental efforts are already pursuing 
the polarization of the CMB; many others will follow soon, 
as the field is literally blossoming. The potential reward from 
this activity is immense, since the CMB is thought to encode 
the solution to several long standing puzzles in Cosmology 
and Fundamental Physics. While the largest contribution to 
CMB polarization (so called E mode) arises due to the ef- 
fect of the scalar perturbations that also seed the large scale 
structure of the Universe, theory predicts that a tiny part 
of the polarized signal is in the form of (yet to be detected) 
B modes. On large angular scales, these modes bear the im- 
print of the stochastic background of gravitational waves 
generated during the inflationary phase of the Universe. At 



the same time, the amount of power in B modes can be 
used to measure the energy scale of inflation, thus probing 
particle physics. Furthermore, the different parity behavior 
of E and B modes opens up the possibility to test for the 
breakdown of fundamental symmetries. 

CM B anisotropies were first discovered by the COBE 
satellite (jSmoot et al.|[l993 ). and highe r resolution experi- 
ment s detected the first a coustic peaks l|de Bernardis et al.l 
I1999I : iHananv et al.l I2OO0I ) in temperature. T he E modes 
were first detected from the ground by DASI (|Kovac et al.l 
l2002h . Soon after, the WMAP satellite produced all sky 
data with a resolu tion down to about 15 arcminutes 
(|Bennett et al.ll2003l ). impacting in particular on the large 
scale E polarization of the CMB, and on the polarized fore- 
ground properties ((Page et al., ,20071 '). The Planck satellit^B 



^ Planck (http://www.esa.int/Planck) is a project of the Euro- 
pean Space Agency - ESA - with instruments provided by two 
scientific Consortia funded by ESA member states (in particu- 
lar the lead countries: France and Italy) with contributions from 



2 Ricciardi, Bonaldi et al. 



l|Tauber et al.|[2009l ) is now measuring the CMB anisotropy 
with an unprecedented accuracy. Lately, experiments are 
focussing on the mapping of small scales total intensity 
anisotropies jReichardt et alj [mOS I . and of polarization 
l|de Bernardis et~ 20091 ') with the ambitious goal of de- 
tecting the B modes for an interesting multipole range. 
The latter projects represent an extraordinary technological 
and scientific challenge, requiring a post-Planck, polariza- 
tion dedicated satellite. 

It must be clearly set forth that building high quality ex- 
periments is not the only necessary condition for the CMB to 
disclose its secrets. It is also of utmost importance to analyze 
the CMB data optimally in order to maximize the informa- 
tion drawn from the data. First, one is looking for a tiny sig- 
nal buried under overwhelming instrumental noise, of both 
statistical and systematic origin. Furthermore, the CMB is 
not the only signal on the sky: other astrophysical sources 
exist, both compact and diffuse, that are powerful emitters 
in the microwave band. Their emission can easily jeopardize 
measurements of the CMB unless a very accurate separa- 
tion of the astrophysical components is achieved. Despite the 
large nu mber of papers on the subject - fromiBrandt et al. 
lll994) tolStompor et all (120091 ) and references therein, see 
iDelabrouille fc Cardosol (|2009l ) for a recent review - the han- 
dles used to achieve component separation are a few. On 
one side, one can linearly combine the multifrequency maps 
in order to extract from the data the most likely compo- 
nent scaling as a bl ackbody (Internal Linear Combination: 
iBennett et al.l [20031 . and references therein), or can exploit 
the statistical independenc e between components (Indepen- 
dent Component Analysis: IStivo iill2006': 'Betoul e et aLll2009l . 
and referen ces therein), or adopt internal template fitting 
procedures l|Martmez-Gonzalez et al"]|2003h . On the other 
side, one can exploit a partial knowledge of foregrounds, 
and in particular of their frequency scalings, in order to pa- 
rameterize them and measure their pa rameters from multi- 
frequency maps, in the pixel dornain (IStompor et al" ] |2009l : 
lEriksen et aLl l2006l: |Bedini et ahlboOSl ) or in the harmonic 



one (IStolvarov et akboOSl ). A crucially important still open 
issue is the estimation of errors on separated maps and their 
propagation through all the steps of the analysis, from the 
determination of spectral properties of the Galactic emis- 
sions to the component maps and CMB power spectra. A 
correct characterization of errors is clearly essential for anal- 
yses of separated maps and for cosmological parameter es- 
timation. In this paper we present a pipeline for component 
separation in polarizatio n based on the Correlated Com- 
pone nt Analysis (CCA, iBedini et al. I l2005l : iBonaldi et al.l 
I2OO6D . including error estimates. Two methods for the esti- 
mation of errors on separated maps are presented and tested 
on realistic simulations of Planck polarization data. 

The outline of the paper is the following. In §[2] we for- 
malize the component separation problem. In §[3]we describe 
our approach: the CCA technique both in the pixel and in 
the harmonic domain, the associated error estimation and 
the reconstruction of the individual components. In §[4] we 
describe the simulations used to test our pipeline. In §[5] 



NASA (USA), and the telescope reflectors provided in a collabo- 
ration between ESA and scientific Consortium led and funded by 
Denmark. 



we assess the goodness of the Galactic spectral indices esti- 
mated on sky patches, used in §|6]to build spatially varying 
spectral index maps. In §[7] we present the component maps 
and the associated error maps. Finally, in §[8]we present and 
discuss the estimated CMB polarization power spectra. 



2 STATEMENT OF THE PROBLEM 

The sky radiation, x, from the direction r at the frequency 
f results from the superposition of signals coming from Nc 
different physical processes Sj : 



(1) 



The signal x is observed through a telescope, whose beam 
pattern can be modeled, at each frequency, as a shift- 
invariant point spread function B{f,v). For each value of 
i/, the telescope defocuses the physical radiation map by 
convolving it with the kernel B. The frequency-dependent 
convolved signal is input to an A'^d-channel measuring instru- 
ment, which integrates the signal over frequency on each of 
its channels and adds noise to its outputs. The output of the 
measurement channel at a generic frequency v is 

Xui-r) = I B{f~f', u') t„{u')sj{f' , v')dr dv' -\-ni,{f), (2) 

where tv{v') is the frequency response of the channel and 
nv{r) is the noise map. The data model in eq. ([2| can be 
simplified by virtue of the following assumptions: 

(i) Each source signal is a separable function of direction 
and frequency: 



(3) 



(ii) B{f,v) is constant within the passband of the mea- 
surement channel containing v. 



These two assumptions lead us to a new data model: 



(4) 



where Bi,{r) is the telescope beam pattern at the effective 
frequency * denotes convolution, and 



Kj = / U{v )fj{u )dv . 



In vector form, we have 



[B*Hs](f) + n(f) 



(5) 



(6) 



where B is a diagonal A^tj-matrix whose elements are Bi,(r), 
H is a constant Nd x Nc matrix whose elements are /i^j, s 
is an Ac-vector whose elements are Sj{r), and n is an A'^^- 
vector whose elements are n^(f). The data model has thus 
become linear and convolutional, with known point spread 
functions. 

Equation ® can be translated to the harmonic domain 
where, for each transformed mode, the linear convolutional 
model becomes linear and instantaneous: 



X = HBS + N, 



(7) 
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where X, S, and N are the transforms of x, s, and n, re- 
spectively, and B is the transform of the matrix B. 

Let us now assume that the beam patterns of the tele- 
scope are the same for all the measurement channels, that 
is, 

B,Ar)^B{r). (8) 
By virtue of this position, eq. Q becomes 

^•^{r) ^^h^iSj{f) +n^{f), (9) 

with Sj{r) = [B * s]{r), or, in vector form, 

X = Hs + n. (10) 

Hereafter we drop the overbar from the symbol of the source 
vector. 

It is worthwhile to note that eq. ((3| comes from an 
important assumption: the spectral properties of the astro- 
physical sources are spatially uniform on the rip pixels con- 
sidered. This assumption must be dealt with carefully be- 
cause the Galactic components are spatially varying, as we 
discuss better in §[4] Our strategy to overcome this difficulty 
is to apply our method separately to sky patches, where the 
foreground properties can be safely assumed to be uniform. 
Generally, the assumption leading to eq. ([8]), needed to build 
an instantaneous data model, is also not realistic. A simple 
way to apply the model [eq. (|10|) ] to a general case is to 
equalize the resolution of the instrumental channels to the 
worse one. For the harmonic domain this is not needed, so 
that the full instrumental resolution of each channel can be 
exploited. 

2.1 Strategy for component separation 

Equations (|10|) and (O clearly show that the key ingredi- 
ent to estimate the source vector is the mixing matrix H. 
The Correlated Component Analysis (CCA), described in 
the next section, gives an estimate of the mixing matrix. 
This estimation could be performed in the pixel domain (§ 
13. 1|) and in the harmonic domain (§[3]2]). Both methods ex- 
ploit a tessellation of the data set and an estimation of the 
mixing matrix patch by patch. Once we have an estimate of 
H, we can compute a suitable matrix W, sometimes called 
reconstruction matrix, allowing us to obtain an estimate s 
of the components from the noisy data x: 

§ = Wx. (11) 

with a reconstruction matrix W = W(H). 

The reconstruction as in eq. (|11|1 could be done in prin- 
ciple both in pixel and in harmonic domain. As we will dis- 
cuss in § I3.4l we choose to perform the reconstruction in pixel 
domain as this technique allows us to account for spatially- 
varying mixing matrix in a more direct way. 



3 CORRELATED COMPONENT ANALYSIS 
(CCA) 

The CCA exploits a second-order statistics to estimate the 
mixing matrix from the statistics of data and noise. To re- 
duce the number of unknowns, the mixing matrix is parame- 



terized through a parameter vector p, such that H — H(p). 
To choose a suitable parameterization for H we use the fact 
that its elements are proportional to the spectra of astro- 
physical sources. As discussed in § U this allows us to re- 
duce substantially the number of unknowns in the mixing 
matrix with respect to the original Nd x Nc elements to be 
estimated. 



3.1 Pixel domain CCA 

From the pixel-space data model in e g. (1101). we easily de rive 
the following second-order statistics (jBedini et akllioosl ): 

Cx(r,V) =HCs(r,V)H^ + Cn(r,V). (12) 

The quantities Cx(t, ?/;), Cs(r, t/;) and Cn{T,tp) are the co- 
variance matrices of data, components and noise, respec- 
tively, computed for a generic two-dimensional shift {T,tp). 
For example, for the data covariance matrix we have: 

C.(r, V) = ([x(e, 77) - m1[x(? + r, r, + V-) - (13) 

where ^ and rj are the coordinates of the two dimensional 
space, T and 1/) are increments in ^ and rj, (...) denotes expec- 
tation under the appropriate joint probability distribution, 
fi is the mean vector and the superscript T means transpo- 
sition. The application to sky patches is straightforward, as 
we simply need to compute the covariances on a suitable list 
of pixels. 

In eq. (|12|) . Cx(t, 1/;) and C„(r, 1/;) can be computed 
from the data and the known statistics of noise, while H 
and Cs(r, ?/)) are unknown. 

Once we consider eq. (|12|) for a sufficient number of 
shift pairs(T, ?/)), both p (and hence H) and Cs(t, V") can be 
estimated by minimizing the functional: 

*[C.,H] = ^|[HC,(r,^)H^-Cx(r,V') + Cn(r,^) 1^(14) 



3.2 Harmonic-domain CCA 

By dropping the equal beam assumption [eq (|8]l] and relying 
on the data model of eq. (O in the harmonic domain, we 
easily derive an equivalent of eg. (|12p in terms of power 
cross-spectra jBedini &: Salerndl2007l '): 

Cx = BHC.H^Bt + C„ (15) 

where B is the transform of the matrix B introduced in 
the previous section, and the dagger superscript^ denotes the 
adjoint matrix. The matrices Cx(^), Cs(/) and Cn(^), all de- 
pending on the multipole £, are the cross-spectra of the data, 
sources and noise, respectively. Formally, they are obtained 
by applying the spherical harmonic transform to the covari- 
ance matrices in eq. (12). When working in small patches, 
the cross-spectra are calculated by averaging circularly the 
2-dimensional discrete Fourier transform on the rectangular 

If we reorder the matrices Cx(i') — Cn(^) and Cs(^) into 
vectors d{£) and c(^), respectively, we can rewrite eq. (15) 
as 

d{£) = ii^{£)c{£), (16) 
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where Hk(^) = [B(^)H] (g) [B(^)H] and the symbol ® de- 
notes the Kronecker product. Hk(^) and c{£) contain the 
unknowns of our problem (namely, the vector p and all the 
components of the source cross-spectra), and d(^) would be 
known if the data cross-sp£ctra were knownjn our case, we 
use binned cross-spectra, Cx(^), Cs(^) and Cn(^), obtained 
by averaging the transforms onto suitable spectral bins Dg. 
Thus, the cross-spectra of data can be estimated from the 
available data samples: 



(17) 



where the pairs are the modes contained in the spectral 
bin denoted by in the transformed domain and Mj; is the 
number of Fourier modes contained in the i'-th spectral bin 
D^, £ — 1, . . . ,£max- The set can be any subset of the 
Fourier plane, such as an annular bin defined by its mean 
radius and its thickness, which can easily be related to a 
specific ^-interval in the spherical harmonic domain. Since 
the left-hand side of eq. (16) can only be evaluated through 
the empirical data cross-spectra of eq. (17), the data model is 
affected by an estimation error e{£), with covariance matrix 
N,(i): 



d(^) = Hk(^)c(^) + e(^), 



(18) 



where d{£) is now computed using the observed data cross- 
spectrum matrix. Vectors e{£) represent the differences be- 
tween the components of the actual data cross-spectrum ma- 
trix and those evaluated from the data through eq. p7p . 
and, of course, are ordered as vectors in the same way as 
Cx(-^). This ordering induces the structure of their covari- 
ance matrices, Ne(i). Let us express the mapping between 
the indices j and k in matrix [Cx(^) — Cn(^)] and the cor- 
responding index i in vector d{£) as follows 



d.(^) = [Cx(^)-C„(£)],(,),fe(,). 



(19) 



With this position, if the estimation errors in d{£) are un- 
correlated to each other, the matrix T^e{£) is diagonal, and 
its entries are given by 



N^Mi) = 24,) /M,-, if j{i) = k{i) 
N.Me) = a|(,)a^(,)/M,-, if j(i) / fc(i) 



(20) 



where is the known variance of the j-th element in the 
noise vector. Details o n this derivation can be found in 
iBedini fc Salerno! (|2008l ). 

Let us now exploit all the significant spectral bins, that 
is, let us assume a suitable set [1, £max] for £, and rewrite eq. 
(^18) by stacking all the quantities for the available values of 
t 



dv = [d(l),d(2),...d(£„ 
cv = [c(l),c(2),...c(^"™, 
€v = [e(l),e(2),...e(W 



HkB 



Hk(l) 








Hk(2) 







Hk(i„ 



(21) 
(22) 
(23) 

(24) 



By these positions, and bearing in mind that matrix HkB 
is completely specified by the parameter vector p, eq. (|18l) 
becomes 



HkB(p) ■ Cv + ev, 



(25) 



which allows us to estimate the parameter vector and the 
source cross-spectra by minimizing the functional: 

*[p,cvl= (26) 
= [dv - HkB(p) ■ cv]^NEB~^[dv - HkB(p) • cv] 
+ XcvCcv 



with 



N,(l) 








N,(2) 







N,(L 



(27) 



The term Ac^Ccv is a quadratic stabilizer for the source 
power cross-spectra whose estimation is an ill-posed prob- 
lem. The matrix C must be suitably chosen and the pa- 
rameter A must be tuned to balance the effects of data fit 
and regularization in the final solution. The functional in 
eq. H26p can be considered as a negative joint log-posterior 
for p and cv , where the first quadratic form represents the 
log-likelihood, and the regularization term can be viewed as 
a log-prior density for the source power cross-spectra. 

3.3 Foregrounds spectral parameters error 
estimation 

A standard, theoretically-based, method to estimate the er- 
rors in the recovery of the spectral parameters, Ap, leading 
to an error in the mixing matrix AH, relies on the anal- 
ysis of the marginal probability for the parameters p. In 
the next section we describe the formalization of this tech- 
nique for harmonic domain CCA. As an alternative to this 
method, we also implemented a simpler technique, exploit- 
ing the redundancy of solutions obtained for different sky 
patches. This allows us to provide an error estimation also 
for pixel domain CCA, for which the formalization of the 
marginal probability method is still under development, and 
to cross-check the results obtained for the harmonic domain 
CCA. 



3.3.1 Marginal probability method for harmonic domain 
CCA 

To evaluate the estimation errors on p, we can first derive, 
from eq. (|26|l . the joint distribution of p and cv and then 
marginalize it by integrating out cv. 

By developing the functional in eq. ([26} and taking its 
exponential, it is easy to see that the marginal density of p 
conditioned to dv is given by 



p(p|dv 



oc / e 



(28) 



-id?;N,„-idv 



which, by dropping unessential constants, becomes 
l|Bedini fc Salernd(2008l ): 

p{p\dv) oc v/det(HkB^N,B"'HkB + AC)-i x (29) 
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gi{N,B"^HkB){HkB^N,B"^HkB+AC)-l{HkB^N,B"^dv) 

Studying the behavior of this marginal distribution is not 
difficult, since p is a low-dimension vector (typically, it has 
two or three components). From eq. (|29p . all the quantities 
related to the parameter distributions can be evaluated and, 
possibly, exploited to estimate all the relevant reconstruction 
errors. As a quantitative index of the estimation errors Ap 
we simply assume the standard deviations evaluated from 
the normalized-exponent marginal posterior. 

3.3.2 Spatial redundancy method 

On a certain patch of the sky, we call p the true spectral 
parameters and p the estimated ones; the quantity we want 
to estimate is Ap = p — p. Let us suppose to have a set 
of parameters estimated on a sample of sky patches, widely 
overlapping to the considered one. For this sample we call 
the true and the estimated spectral indices {pj} and {pj} 
respectively. 

For the moment we make a simplifying assumption, 
which we will relax later: the spectral behavior of the com- 
ponents is uniform over the area covered by the sample of 
patches we are considering. In this case we have: 

{P.} = P, (30) 

and the estimated parameters {pj} are different measure- 
ments of p. Thus, we can estimate Ap by comparing the 
actual estimation in the considered patch p with the expec- 
tation value of p computed on the sample of patches: 

Ap=({p.})-P. (31) 

In a realistic situation, however, the spectral indices are 
spatially-varying on scales smaller than that of the patches. 
Thus, for each patch the true spectral indices have a certain 
distribution. Our assumption is that in this case eq. pop . 
though not strictly true, still approximately holds. In fact, 
due to the autocovariance of the foreground signal, we do not 
expect discontinuities in the spectral index distribution in 
nearby regions. Moreover, as the patches are partially over- 
lapping, their mean spectral indices are highly correlated. 
Thus, we assume that eq. (|3H) can still be used for an ap- 
proximate error estimation, provided that the patches are 
small and overlapped enough. A test of this method using 
simulations is described in §(5] 

3.4 Reconstruction of the components 

Adopting the linear mixture model in eq. (|10p . we can find 
a solution of the component separation problem of the form 
of eq. pip . By exploiting the CCA results on sky patches 
we are able to synthesize spatially varying spectral index 
maps (see §[6] for details). Once we have an estimation of 
the mixing matrix H, a suitable choice for the reconstruction 
matrix W is the Generalized Least Square solution (GLS): 

W = [H'^C-^H]-'h'^C„\ (32) 

which only depends on the mixing matrix and on the noise 
covariance Cn- By applying eq. (|32p in pixel space we are 
able to preserve the spatial variability of the estimated mix- 
ing matrix. This local information is very valuable, as the 
spectral properties of the foregrounds depend on the line 



of sight. One disadvantage of using the pixel domain is that 
the noise covariance Cn among different pixels does not van- 
ish, so that the full calculation of eqs. H32p and pil) is very 
computationally demanding. For full resolution maps, hav- 
ing ~ 10^ pixels, computing the full Cn is infeasible in prac- 
tice, and we are forced to take into account only the diagonal 
noise covariance, i.e. to neglect any correlation between noise 
in different pixels. The full calculation can be performed on 
low resolution maps having ~ 10'' pixels. Such code, un- 
der development, will provide a full noise covariance for the 
reconstructed low resolution CMB map and will naturally 
couple CCA with low resolution power spectrum estimators 
such as Bol-Pol l|Gruppuso et al.l 12009!). Another disadvan- 
tage of the pixel domain is that to apply the simple relation, 
eq. (|32p . we need to equalize the beams, thus losing part of 
the instrumental resolution. This can be avoided again at 
the cost of an increased computational complexity. An it- 
erative multi-resolution solver is under development, which 
will play perfectly with the harmonic CCA fully exploiting 
the native channel resolution. 



3.5 Errors on component maps 

In general, each component map estimated through eq. (|lip 
is affected by both the instrumental noise and the residual 
contamination from the other components. The former has 
covariance matrix 

cov(s) = WCnW" (33) 

The latter will be estimated propagating the errors on spec- 
tral parameters to the actual maps of individual compo- 
nents. This is done by performing separation using different 
realizations of the spectral parameters from their associated 
distribution characterized by the uncertainties described in 
this Section. Following this idea, in §[7] we present an ap- 
proximated analysis and we demonstrate that, in our sim- 
ulations, the error propagation from foregrounds spectral 
parameters to component maps is satisfying. Notice that a 
straightforward, although computationally expensive, exten- 
sion of the present treatment could be exploited to propagate 
covariances of the recovered errors on spectral parameters, 
by simply computing them while conducting spectral indices 
estimation in parallel for all patches. This step might be cru- 
cial especially at large scales, where foregrounds do have cor- 
relations which could impact significantly on the separation 
errors in the form covariances between the distribution of 
spectral parameters associated to nearby patches. This im- 
plementation is on the other hand beyond the scope of the 
present paper that regards small and intermediate scales but 
we plan to implement it for the pipeline of real data analy- 
sis. 



4 SIMULATED SKY 

The main diffuse components present in the Planck chan- 
nels are the CMB and the emissions due to our own Galaxy, 
namely thermal and anomalous dust, synchrotron, and free- 
free. Free-free and anomalous dust emission are expected to 
be essentially unpolarized. Since we deal here with polariza- 
tion data, we are left with CMB, synchrotron and thermal 
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Table 1. Central frequencies, u, and angular resolutions for the 





Figure 1. Spectral index maps of synchrotron (top) and dust 
(bottom) components, used in the simulated sky 



dust emissions only. We note, however, that the method de- 
scribed in this paper can be straightforwardly applied also 
to total intensity data. 

Our simulation of the diffuse emission exploits all the 
available information from existing public pre-Planck sur- 
veys. The simulated CMB map is based on a standard 
ACDM model consistent with WMAP 5-years cosmological 
parameters, with a tensor to scalar ratio r = 0.1. The rms 
fluctuations of the CMB are expressed in antenna tempera- 
ture, rA,cMB(j'), whose frequency dependence writes: 



^ , , {hu/kTcMs)^ exp{hu/kTcMB] 

(exp(fti//fcTcMB) - 1)^ 



(34) 



where h is the Planck constant, k the Boltzmann constant, 
TcMB = 2.726 K. 

The simulation of the diffuse Galactic emissions is based 
on lMiville-Deschenes fc Boulangeil (|2007l ) and implemented 
by the Planck working group on Component Separation. 
The starting point towards building the synchrotron emis- 
sion template is an extrapolation of the 408 MHz map of 
iHaslam et al. I l| 19821 ) from which an estimate of the free- 
free emission has been removed. Due to the poor resolution 
of the Haslam map (52 nominal arcmin) small structures 
have been artificially added using the p rocedure presented in 
iMiville-Deschenes fc Boulangeil l|2007l ). The fraction of po- 
larization is derived from a model of the magnetic field in- 
cluding spiral and turbulent components based on WMAP 
5yr results. On large scales (> 5deg) the polarisation an- 
gle is WMAP constrained; on smaller scales it relies on the 
Galactic magnetic field model. The spectrum, in antenna 
temperature, is assumed to follow a power law: 



Planck channels considered in the present 


study 






u(GUz) 30 44 70 100 


143 


217 


353 


FWHM (arcmin) 33 24 14 9.5 


7.1 


5.0 


5.0 



rA,synch(j^) OC U 



(35) 



where the synchrotron spectral index jSs varies with the po- 
sition on the s ky. To this aim, wc use the map of spectral in- 
dices given bv lGiardino~e t al. (2002'). This map shows struc- 
tures on scales of up to 10 degrees, with /3s varying from 2.5 
to 3.2 (see the top panel of Fig. [l|. 

The simulation of the polarized therma l dust emission is 
based on model 3 of lFinkbeiner et "all (|l999l). This model ex - 
trapolates the 100 /^m brightness map ( Schlegel et al.|[l998l ) 
assuming grey body spectra: 



rA,dust(l^) oc 



./9d + l 



ip{hiy/kTiust) - 1 ' 



(36) 



with Tdust = 18 K and spatially varying emissivity index fid, 
obtained fro m the map of 1 00/24 fim flux density ratios 
published by ISchlegel et all (|l998l '). The map of Pd. shows 
structure down to arcminute scales (see the bottom panel of 
Fig.[T} with Pd varying from 1.44 to 1.6. The polarized inten- 
sity is obtained multiplying the brightness map by a polar- 
ization fraction map extracted from WMAP 5yr data with 
the help of the model of the previously mentioned Galac- 
tic magnetic field model. In this simulation the polarization 
fraction is not frequency dependent. 

Component maps have been produced at central fre- 
quencies of all Planck channels. They are then co-added 
and smoothed with nominal Gaussian beams (see Table 
[l}. A white inhomogeneous noise is synthesized using the 
block diagonal part of the predicted noise covariance matrix 
given the Planck nominal integration time (14 months) 
and scanning strategy. The correlation between noise in dif- 
ferent Stokes parameters for each pixel has also been re- 
produced. This dataset has been complemented with the 
simulated WMAP 23 GHz map after 5 years of survey. This 
map has a resolution of 58 arcmin and the noise is simulated 
starting from the diagonal noise covariance provided by the 
WMAP team (http://lambda.gsfc.nasa.gov). This ancillary 
product is used to help tracing the low frequency foreground 
component for the CCA run, as described in §[3 



5 ESTIMATING THE MIXING MATRIX 

PARAMETERS AND ERRORS WITH CCA 

We applied both pixel and harmonic domain CCA to the 
polarized Planck simulations described in the previous sec- 
tion (7 frequencies from 30 GHz to 353 GHz), with and with- 
out the ancillary 23 GHz WMAP channel. For the harmonic 
domain version we use the channel maps at their native reso- 
lution. In the application to polarized data we perform sep- 
arated runs of the Q and U maps so the c(^) in eq. ([18} 
simply represent cross-power spectra. We can evaluate the 
local mixing matrix using both information from Q and U 
maps. In practice, due to the lower foregrounds signal in 
the U maps, the mixing matrix is basically defined by the Q 
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maps analysis only. The application of pixel-domain CCA 
requires that all the channel maps have the same resolution. 
Thus we preliminarily degraded the maps to the resolution 
of the lowest frequency channel (33 arcmin to work with 
Planck data only, 58 arcmin to include WMAP 23 GHz 
channel) convolving the maps with a Gaussian beam. The 
data model includes three diffuse components: CMB, syn- 
chrotron and thermal dust, parameterized as in eqs. (|34|) . 
(|35p and psp respectively. We estimated two free parame- 
ters, the synchrotron and the dust indices, which were al- 
lowed to vary in the ranges 2.3 ^ /3s ^ 3.5 and 1. ^ Pd ^ 2.5, 
while we kept Tdust ~ 18 K. We note that we do not explic- 
itly account for any modeling error in our analysis, as the 
mixing matrix parametrization assumed by CCA exactly re- 
flects the one exploited for the data generation. Thus, the 
only systematic error is given by the assumption that the 
spectral properties are constant within sky patches. Includ- 
ing the effect of an imperfect modeling is beyond the scope 
of this paper, whose goal is to evaluat e CCA performances 
per se. The problem was addressed in lBonaldi et al.l (|2007l ) 
where several models for the anomalous emission were tested 
for the analysis of WMAP temperature data with CCA. In 
any case, we believe that such model uncertainties should 
be less severe in polarization. 



The choice of the patch size for the CCA run is a trade- 
off between the need to have uniform foreground properties, 
which calls for small patches, and the need to have enough 
statistics, which calls for bigger ones. The latter is obvi- 
ously related to the instrumental resolution of the data, as 
the statistics is ultimately determined by the number of res- 
olution elements in the considered region of the sky. In the 
case of harmonic-domain CCA we have the additional con- 
straint of the maximum patch size allowed by the planar 
approximation. However, the possibility of this version of 
CCA to handle frequency-dependent instrumental beams, 
and thus to exploit the full resolution, generally allows the 
use of smaller patches compared to pixel domain CCA. In 
this work we divide the sky in square patches for both pixel 
domain and harmonic domain CCA. The pixel based version 
does not suffer of any constraint regarding the shape of the 
patch. The current harmonic version instead can work only 
on square regions. We adopt a patch size of 40° x 40° for 
pixel-domain CCA, of 30° x 30° for harmonic-domain CCA. 
The centers of the patches are equally spaced in latitude and 
longitude with shifts of 3°, up to a maximum central lati- 
tude of ±30° so that we have 2520 patches in the sky. This 
helps us to build a smooth spectral index map and allows 
a more localized spectral index estimation, as described in 
§[6l Moreover this provides the redundancy needed for error 
estimation as described in ij |3.3.2l for which we considered 
samples of patches overlapping by more than 60%. We note 
however that this purely geometrical partition of the sky is 
not driven by any astrophysical reason. If we could achieve 
a partition that maximizes the uniformity of the spectral 
properties within each patch, the CCA estimation would be 
more accurate. The full run took ~ 10 hours on a parallel 
machine at NERSC using 120 processors. Convergence has 
been reached in ~ 80% of the patches for pixel domain CCA 
and ~ 90% of the patches for harmonic domain CCA. 



5.1 Evaluation of the results 

5.1.1 Actual errors on parameter estimation 

Ideally, to evaluate the quality of CCA results we should 
compare the true synchrotron and dust spectral indices with 
the estimated ones for each position in the sky. However, 
CCA only provides spectral indices per patch, and the true 
spectral indices in general vary within the patch. Thus we 
computed: 

Ap = p - p, (37) 

where p = [Ps, Pd\ are the flux weighted averages of the true 
parameters over the patch and p = [Ps, Pd\ are the esti- 
mates provided by CCA. This choice for the "true" spectral 
indices per patch seems to be the most appropriate because 
the CCA results mostly reflect the spectral properties of 
the brightest foreground structures in the considered patch. 
Such bright structures are obviously those that need to be 
more accurately removed to produce a clean CMB map. In 
any case, a systematic difference is expected because the 
"true" mean spectral indices and the ones recovered by CCA 
have somewhat different meanings. 

In Fig. [2] we show for each parameter and for both 
pixel domain and harmonic domain CCA the distribution 
of "true" errors for the full sample of patches obtained from 
simulated Planck Q maps with and without the simulated 
WMAP 23 GHz channel. Table [2] reports the offset (from 
zero) of the mean and the standard deviation of a gaussian 
fit computed for each histogram. We note that the offset is 
in some cases comparable with the rms; this systematic error 
is induced by the spatial variability of the true indices, as 
mentioned above. We verified that, once we adopt spatially- 
uniform spectral indices for the data generation (/3s = 2.9 
and Pd = 1-7), the offset disappears while the width of the 
distributions is almost unchanged. Errors in the recovery of 
the synchrotron spectral index are bigger than those for the 
dust. This is due to the fact that the dust component is 
much better traced by the frequency coverage of Planck. 
The inclusion of the WMAP 23 GHz channel almost cancels 
the offset from zero of the mean error on the synchrotron 
spectral index for harmonic domain CCA, while slightly in- 
creasing the rms value (see Table [2|. For pixel domain CCA 
the advantage of a broader frequency range yielded by the 
inclusion of the 23 GHz map is more than compensated by 
the degradation of the resolution of the whole data set; as a 
result, both the mean offset and the rms value of errors some- 
what increase. In the following we only consider the results 
obtained with Planck maps alone for pixel domain CCA, 
and with Planck -|- 23 GHz data for harmonic domain CCA. 
The rms errors in the synchrotron and dust spectral indices 
are 0.08 and 0.019 respectively for pixel domain CCA, 0.045 
and 0.009 respectively for harmonic domain CCA. Thus the 
harmonic domain CCA performs better than the pixel do- 
main version. 

5.1.2 Actual versus estimated errors 

As illustrated by Fig.[3]the errors estimated via the marginal 
distribution method (§[3]3Tl]) are well correlated to the 
"true" errors above A/?s — 0.05 or A/?d — 0.01. When the es- 
timated errors are small, the "true" errors are small as well, 
although the two quantities are not correlated. This is due 
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Figure 2. Distribution of errors in the recovery of synchrotron (left) and dust (right) spectral indices for the full sample of sky patches 
using only the simulated Planck maps (grey) and adding the simulated WMAP 23 GHz map (black). Top: pixel domain CCA; bottom: 
harmonic domain CCA 




Actual error Actual error 



Figure 3. Absolute values of estimated vs actual error on synchrotron (left) and dust spectral index (right) for the harmonic domain 
CCA; error estimation through the marginal distribution method. Dots are the values for individual patches, squares and error bars are 
the median values and quartiles computed for groups of 100 points, respectively 



Table 2. Offsets of the mean and rms values for the error distri- 
butions in Fig. [2] (Planck dataset/Planck + 23 GHz data set) 







offset 


rms 


Synchr. 


Pixel CCA 


-0.009/0.016 


0.080/0.107 




Harmonic CCA 


-0.027/0.002 


0.036/0.045 


Dust 


Pixel CCA 


-0.005/-0.004 


0.019/0.022 




Harmonic CCA 


0.005/0.003 


0.010/0.009 



to the effect of noise which hampers the calculation of the 
marginal probability distribution with sufficient accuracy. 
The largest errors are moderately underestimated; this is 
not particularly worrisome however, since large errors corre- 
spond to regions where the foreground signals are very weak 
and can be removed with sufficient accuracy adopting the 
mean parameter values determined in the rest of the sky. 
The spatial redundancy error estimation method (S I3.3.2|I 
applied to both pixel and harmonic CCA yields estimated 
errors systematically lower than the "true" ones. However, 
the latter are well matched by the third quartiles of the 
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Figure 4. Rms errors in the estimation of spectral indices as 
a function of Galactic latitude. Solid lines: full sample; dashed 
lines: final sample obtained with spatial redundancy error esti- 
mation method; dotted lines: final sample obtained with marginal 
distribution error estimation method. The errors increase at high 
Galactic latitudes, where the foreground signal are weak and there 
removal is therefore not a critical issue. 



distribution of estimated errors. As expected, the spatial re- 
dundancy method performs better for the harmonic than for 
the pixel domain CCA, thanks to the higher spatial resolu- 
tion (30° X 30° versus 40° x 40° patches). 

Despite the non-idealities mentioned above, errors es- 
timates produced by both methods are good enough to be 
exploited in the next step of our pipeline. In the following 
analysis we will restrict ourselves to values of each spectral 
index j3 within the range |A/3| < /i^^ + cr^^, where A/3 is 
the estimated error on j3 and and are the mean 
and the standard deviation of the estimated errors over all 
the patches. This set of values constitutes the final sample. 
Figure |4] compares the rms values of the actual estimation 
errors as a function of Galactic latitude for both the full 
and the final sample, showing that indeed this procedure 
re moves the most discre pant values. As already pointed out 
bv lBonaldi et al] (|2007l ). the performance of CCA is highly 
dependent on the intensity of the foreground signal in the 
considered patch. Estimations of the spectral parameters are 
better close to the Galactic plane and worsen with increas- 
ing Galactic latitude where, however, the foreground signals 
are anyway very weak and their removal is therefore not a 
critical issue. Our method allows us to flag the worse esti- 
mates and to replace them with average values determined 
in the rest of the sky. 



6 BUILDING FULL-SKY MAPS OF 

SPATIALLY- VARYING SPECTRAL INDICES 

The information on spectral parameters provided by CCA 
for sky patches can be combined to build full-sky maps 
of the estimated, spatially-varying spectral indices, and 
the associated estim ation error maps. We use the Healpix 
ijGorski et al.l l2005l i pixelization, with resolution parame- 
ter NSIDE=64, corresponding to pixel areas of about 1 
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Figure 5. Dependence of the weight wl on the distance from the 
patch center normalized to half of the patch size A 



square degree. Let us consider the component p (e.g. the 
synchrotron spectral index) of the parameter vector p. Due 
to the partial overlap of sky patches, the j-th pixel of the 
Healpix map, whose position on the sky is centered at rj, 
belongs to several different patches for each of which we 
have obtained an estimate of p. We can therefore compute a 
weighted mean of the values of the parameter pj and of the 
error Apj as follows: 



Pj = 



Apj = 



p(i-0 ■ wl(ri, rj) ■ w2{ri 
J2'iLi^Hri,rj) ■w2{ri) 



(38) 
(39) 



where the sum is over all patches in the final sample, defined 
in 5 15. 1.21 and wl and w2 are weight functions. The former, 
'wl{ri,rj), depends on the distance between the j-th pixel 
and the center of the i-th patch, ||ri — rj||, normalized to 
the half size of the patch, A/2. We used the function shown 
in Fig. [5] equal to 1 for | |ri — 1 1 ^ A /2, and equal to for 
||ri - TjII > A/2. The weight io2(ri) = ii)2(Ap(ri)) depends 
on the estimated error for the i-th patch trough the relation: 

Ap(ri) 



!;2(r0 = 1 - 



<Ap(rO)' 



(40) 



approaching 1 as Ap(ri) — >■ and approaching when 
Ap(ri) — >■ max(Ap(ri)), max(Ap(ri)) being the maximum 
error estimated for the final sample. 

We can associate to the j-th pixel of the Healpix map, 
whose position on the sky is rj, the parameter pj and the 
error Apj as follows: 

Our spectral index and error maps are undefined in re- 
gions, mostly at high Galactic latitudes (|6| > 30°), where 
the Galactic emissions are too low and the CCA cannot 
produce reliable estimates of their spectral parameters. For 
these regions we have adopted the mean values of the pa- 
rameters found for the rest of the sky, using a smoothing 
function to soften edge effects. 
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Figure 6. Distribution of residual (true-estimated) synchrotron (left) and dust (right) spectral index maps in the latitude range —30° < 
b < 30°. The pixel size used is ~ 0.84 square degree (nside 64). Top: pixel domain CCA; bottom: harmonic domain CCA. 



6.1 Evaluation of the quality of the results 

To evaluate the quality of the maps of estimated spectral pa- 
rameters, we looked at the residual spectral index maps. To 
compute those maps, we regridded the true input spectral 
index maps to NSIDE=64, the resolution of the estimated 
maps. In Fig. [6] we show the distributions of residuals (true 
minus estimated spectral indices) with a Gaussian fit su- 
perimposed. The histograms only contain the pixels in the 
latitude range —30° < 6 < 30° where the spectral indices 
have been estimated. Overall, the results are encouraging. 
Again the harmonic domain CCA works better than the 
pixel domain version for both dust and synchrotron spectral 
parameters. Both approaches give negligibly small mean off- 
sets from the true values, but the harmonic CCA has smaller 
dispersions of the residuals [gs — 0.08, (Jd — 0.01) than the 
pixel domain CCA [gs ~ 0.11, ~ 0.02). 



7 RECONSTRUCTION OF THE 
COMPONENTS WITH GLS 

Having shown that the harmonic domain CCA is superior 
to the pixel domain version, we have used its estimates of 
synchrotron and dust spectral indices to build the spatially- 
varying mixing matrix H and the reconstruction matrix W 
given by eq. (|32|l . To recover each emission component we 
have applied the reconstruction matrix [see eq. (|lip ] to the 
simulated data for Planck channels in the frequency range 
from 70 to 217 GHz, where the CMB is the least contam- 
inated by foreground emissions. On the other hand, just 
for the same reason, this frequency range is not optimal 



for the recovery of foreground emissions. All the maps were 
smoothed to the resolution of the 70 GHz channel (14 ar- 
cmin FWHM). We also performed a run at a 60 arcmin res- 
olution to reduce the contribution of instrumental noise, so 
that we can better appreciate the effect of component sep- 
aration errors on the reconstructed maps. Besides the com- 
ponent maps [eq. (|lip ]. we also computed the variance due 
to instrumental noise [eq. (|33|) ] and to residual foreground 
contamination. The latter is estimated by propagating the 
errors on the mixing matrix parameters to the separated 
components as descibed in § 13.51 In particular we assume 
that the spatial correlation of separation errors in the spec- 
tral index in each pixel are the same as those in the spectral 
index map itself. Thus, in order to evaluate and propagate 
the error in the outputs, we construct Monte Carlo varia- 
tions of the estimated spectral index maps. In doing this, 
we further simplify our treatment by only taking into ac- 
count the two-point correlation function of the map. 

Therefore, if f5s and Pd are the maps of the esti- 
mated synchrotron and dust spectral indices respectively, 
and err Pa, err fid the corresponding estimated error maps, 
for the i-th run we obtain the perturbed spectral index maps 
PI and Pd as follows: 

/3i = /3, + A« • err/3, (41) 
/3i = /3d + Ad-err/3d, 

where As and are two maps with zero mean and unit rms 
synthesized from the 2-point correlation function extracted 
from the corresponding estimated spectral index map. 

For each set of fake spectral index maps obtained 
through eqs. (I41|) we performed the source reconstruction 
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Figure 7. Correlations between the true and the reconstructed Q and U maps of Galactic components for ~ 1° pixels. Upper left 
panel: dust Q (correlation coefficient ~ 0.99); upper right panel: synchrotron Q (correlation coefficient ~ 0.85); lower left panel: dust U 
(correlation coefficient ~ 0.98 ). The synchrotron U (lower right panel) was not detectable. 
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Figure 8. True rms fluctuations of input polarized components 
at the reference frequency of 70 GHz as a function of Galactic 
latitude for a 60' resolution (solid line) compared with estimated 
error due to instrumental noise (dashed line) and estimated er- 
ror due to residual contamination for the marginal distribution 
(dotted line) and spatial redundancy (dashed-dotted line) error 
estimation methods. Upper panels: Q Stokes parameter; lower 
panels: U Stokes parameter. 



by GLS, thus obtaining as many sets of output components. 
We computed the variance due to separation for a certain 
component as the variance of all the results for that compo- 
nent. 



In Fig. [8] we show for all the components reconstructed 
with 60' resolution the rms of the true input component 
as a function of Galactic latitude, compared to the a due to 
noise and imperfect separation, computed from the variances 
output by our method. Even at this low resolution instru- 
mental noise dominates except close to the Galactic plane 
where the S/N ratio is ~ 10. The error estimation yielded 
by the marginal distribution method is more conservative 
than the one obtained with spatial redundancy method; the 
two estimates typically differ by 30%. 

7.1 Quality of the reconstructed component maps 

In Fig. [7] we compare the recovered Q and U maps of 
synchrotron and dust emission with NSIDE=64 (resolution 
~ 1°). The agreement is almost perfect for dust (correlation 
coefficient of 0.98-0.99), and very good for the synchrotron 
Q (correlation coefficient of 0.85). The synchrotron U map 
could not be reconstructed because of the very low S/N ratio 
(see Fig.[Sl) due to the choice of frequency channels, that left 
out the low frequency ones (30 and 44 GHz, and WMAP 23 
GHz) where the synchrotron is much stronger. 

As another figure of merit, we computed the normalized 
correlation multipole by multipole between the true and re- 
constructed map: 



(Cf"'=)(Cf™>' ^ ^ 

where C|'^"'^ and C|*™ are the auto-spectra of the true and 
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Figure 10. True and reconstructed maps for dust (top) and synchrotron emission (bottom). 



the reconstructed component, and C"°''° the cross-spectrum 
between the two. Figure [9] shows that, for dust, the cross- 
correlation is close to unity on all the relevant scales. The 
same is true for synchrotron, on scales larger than a few de- 
grees; on smaller scales the correlation drops because the in- 
strumental noise dominates. In Fig. [10] we show the compar- 
ison between the true and the reconstructed Q synchrotron 
and dust emission maps. 

7.2 Assessment of the component separation 
error estimations 

As outlined above, our method outputs two rms error maps, 
one for instrumental noise and the other for component sep- 
aration errors. The actual error map for each component is 
simply computed as the residual map, i.e. the difference be- 
tween the reconstructed map of a component and the true 
one at the same resolution. This residual map contains both 



noise and component separation errors, and thus has to be 
compared to the global estimated error map, obtained sum- 
ming the variances of the two estimated contributions. 

In Fig. [11] we show the standard deviations of the resid- 
uals as function of Galactic latitude (diamonds) compared to 
the total error yielded by both our error estimation methods 
for the 60 arcmin resolution (lines). The results are gener- 
ally very good, as the total error is correctly predicted. On 
the other hand, since we are generally noise dominated, this 
figure cannot tell much about the goodness of our compo- 
nent separation error estimates, except for dust close to the 
Galactic plane, where the error estimate turns out to be very 
successful. We note however that, as we are using a linear 
estimator [eq. (|lip ] to reconstruct the components, our fi- 
nal results can be decomposed into a noiseless term, only 
affected by component separation errors, plus a noise term: 

s = Wy = W[Hs] + Wn. (43) 
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Figure 11. Comparison of true (diamonds) and estimated rms errors as a function of Galactic latitude at 70 GHz for Q (upper panels) 
and U (lower panels) maps. Dotted lines: marginal distribution method (ii l3.3.1ll : dot-dashed lines: spatial redundancy method (ii l3.3.2ll . 
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Figure 12. Noiseless case: comparison of true (diamonds) and estimated rms errors as a function of Galactic latitude at 70 GHz for 
Q (upper panels) and U (lower panels) maps. Dotted lines: marginal distribution method (^ 13.3.11 1: three dots-dashed lines: spatial 
redundancy method fii B. 3.21 1. At high Galactic latitudes, where both our methods overestimate the component separation errors, such 
errors are irrelevant compared to uncertainties due to instrumental noise. 
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Figure 13. Full-sky Q (left) and U (right) CMB maps obtained by CCA for the noiseless case. 



Here Hs is a noiseless dataset, obtained exactly as described 
in §|4] but without adding noise. By combining the noise- 
less dataset with the same reconstruction matrix estimated 
in the noisy case, we obtain a set of reconstructed compo- 
nents having the same component separation errors as be- 
fore, but without noise. This allows us to test the quality of 
component separation error estimates. The results, shown 
in Fig. 1121 are very encouraging. Even if component separa- 
tion errors are in general highly subdominant, the marginal 
distribution error estimation method is able to correctly es- 
timate them at low and intermediate Galactic latitudes (at 
high Galactic latitudes the component separation errors are 
overestimated but are anyway irrelevant compared to errors 
due to noise) . On the other hand the spatial redundancy er- 
ror estimation method is occasionally underestimating the 
true errors at low latitude. In Fig. [13] we show, as an exam- 
ple, the CMB Q and U reconstructed maps in the noiseless 
case. A visual inspection does not reveal the presence of 
residual Galactic contamination except for a tiny strip on 
the Galactic plane. 



8 CMB POWER SPECTRUM ESTIMATION 

To assess the quality of the Stokes Q and U CMB maps 
obtained with the CCA component separation we compare 
their estimated angular power spectra (APS) to the input 
model used to generate the simulation. In doing so we prop- 
agate to the power spectra the component separation errors 
described above. We employ the ROMAster co de, a pseudo- 
Ci e stimator based on MASTER approach dHivon et al 



2002 ) and extended to cross-pow er spectra llPolenta et al 
2005 ) and polarization (see e.g. iKogut et al.l |2003| . for a 



similar formalism). It is well known that the pseudo Ci 
approach to the CMB power spectrum estimation is sub- 
optimal for the lowest mult ipoles where other tech niques are 
more appropriate (see, e.g.. lGruppuso et al.ll2009l '). However, 
a pseudo-Cf estimator is enough for our purpose of assessing 
the quality of the reconstructed CMB polarization maps in 
the presence of noise and component separation errors. 

We exclude from the analysis the regions that are most 
contaminated by residual foreground contributions as esti- 
mated in the previous section. For this purpose we build a 
mask based on our reconstruction errors, flagging all pixels 
where the sum of the variance errors on the CMB Q and 



U maps is greater than the mean value of the same quan- 
tity across the whole map. The resulting mask is shown in 
Fig. [14] and excludes less than 10% of the sky. 

Having only one final map per astrophysical compo- 
nent, we do not rely here on a cross spectrum analysis but 
on an auto-spectrum approach. This is rather general, and 
achieves a lower fina l noise variance than the cross spectrum 
approach (see, e.g.. [Polenta et al]|2005h . The drawback is 
that we need to model and subtract a noise bias in the data. 
To this extent, we computed the noise bias on the CMB EE 
power spectrum by means of 1000 simulated noise maps. To 
obtain each of them, we simulated one noise map for each 
channel included in the reconstruction of the CMB (70, 100, 
143 and 217 GHz), equalized the resolution of all channels to 
14', and combined them with the reconstruction matrix W 
as described in the previous section. ROMAster uses these 
Monte Carlo data to subtract the noise bias, as well as to 
estimate errors on the APS due to noise by computing the 
empirical variance of the realization. 

To compute the error bars due to residual foreground 
contamination, we produced a further set of 100 CMB maps 
by perturbing the input spectral index maps as described 
in the previous section. For each of them we repeated the 
computation of the power spectrum and corrected for the 
noise bias, relying for the latter purpose on a smaller (~ 10) 
set of noise-only maps. The noise bias has been estimated 
each time with the reconstruction matrix used to obtain 
the corresponding CMB map. Even if quite computation- 
ally demanding, this procedure is needed because we are in 
the noise-dominated regime, and a small error in the noise 
bias can substantially affect the estimation of the polarized 
power spectrum. Once we got our 100 unbiased CMB power 
spectra, we finally computed the errors due to component 
separation as the standard deviation of the sample for each 
considered multipole bin. 

In Fig. [15] we show the noise and the component sep- 
aration error bars compared to the EE power spectrum of 
the fiducial model. As we can see, the noise contribution is 
dominating even on the smallest multipoles over the com- 
ponent separation error, which is at most a small correction 
to the error budget. In Fig. [16] we show the EE power spec- 
trum estimated in a realistic Planck case (diamonds); the 
la errors are shown by the shaded area. The results for the 
noiseless case (squares with the component separation error 
bars) show that the accuracy of our estimation of the power 
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Figure 14. Mask used for the computation of the EE power 
spectrum, excluding 9% of the sky on the basis of the estimated 
variance due to component separation on the CMB Q and U 
maps. 




Figure 15. EE power spectrum of the input model (solid) com- 
pared to the rms errors due to noise (dot-dashed line), separation 
(dashed line) and total (solid line) . The dotted line represents the 
associated cosmic variance error. 



spectrum is not limited by component separation but rather 
by the effect of the noise. The overall quality of the recovered 
spectrum is impressive, especially when compared to the to- 
tal Galactic emission at 70 GHz, shown on top of the plot. 
Although a discussion on the detectabil ity of CMB B-modes 
is beyond the scope of this paper (see iBetoule et al.l l2009l : 
lEfstathiou fc Gratton|[2009l . for recent analyses), we briefly 
address here the potential of the GCA in this res pect. An 
interesting point, raised by lEfstathiou et al.l (|2009l ). is that 
some foreground subtraction methods, such as the ILC, are 
not well-suited for applications to B-mode detection, as they 
introduce an offset caused by the dominant E-mode polar- 
ization pattern which prevents estimating B modes even in 
absence of noise. To check whether this also applies to our 
component separation pipeline, we estimated the CMB B- 
mode power spectrum from the CMB maps recovered by 
the CCA in the noiseless case. It is important to note, how- 
ever, that this cannot be really considered as the application 
of our pipeline to an ideal (noiseless) experiment because, 
while the component reconstruction has been made by set- 
ting to zero the noise term in eq. (|43p . the spectral param- 
eters going into the reconstruction matrix used to perform 
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Figure 16. Recovered EE power spectrum in the realistic Planck 
case (diamond) and 1 a uncertainties (shaded area) ; recovered EE 
power spectrum in the noiseless case (squares) with 1 a compo- 
nent separation error bars; EE input CMB map (gray solid line). 
On top of the plot we also show, for comparison, the total Galac- 
tic emission of our simulation at 70 GHz (solid line); the dashed 
line includes the effect of the 70 GHz Planck noise. 



component separation are those estimated by assuming the 
real pre-flight estimates of Planck noise level for the nomi- 
nal mission duration (14 months). The errors on foreground 
spectral parameters, and therefore those on the recovered 
CMB map, are thus substantially larger than those expected 
for an ideal noiseless experiment. 

The results are shown in Fig. 1171 For the tensor to scalar 
ratio adopted in the simulations (r = 0.1), our component 
separation method proves to be capable of recovering the 
B-mode power spectrum over a quite broad multipole range 
without any noticeable spurious feature due to residual fore- 
ground contamination or to cross-correlation with the E- 
mode. Note that the error bars due to component separation 
shown in Fig. [T7]are estimated in the presence of noise while 
the CMB recovery is made using noiseless maps. Since the 
noise dominance will be much stronger for the B-mode than 
for the E-mode, the subtraction of the noise bias will be a 
very delicate and computationally demanding, but concep- 
tually easy, operation. 



9 CONCLUSIONS 

We presented and tested on realistically simulated Planck 
polarization data a p ipeline for component separation based 
on th e CCA method (|Bedini et al.ll2005l : lBonaldi et al.ll2006l . 
l2007f) . This method exploits the data statistics to estimate 
the frequency behaviour of the diffuse components superim- 
posed to the CMB, described in terms of a limited set of 
parameters, which in our case are synchrotron and thermal 
dust spectral indices. 

The most recent harmonic domain version of CCA 
proved to be even superior to the original pixel domain one, 
whic h however gave very good resu lts for simulated Planck 



(Bonaldi et al 



(jBonaldi et al 



.20061 : iLeach et al] .2008 ) and WMAP data 
20071 ') temperature data. As another improve- 



ment, we worked out a method to combine the results ob- 
tained with CCA on several regions of the sky to provide 
spatially-varying maps of spectral indices, with rms errors 
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Figure 17. Recovered B-modc power spectrum in the noiseless 
case (squares with la component separation error bars) compared 
with the input model (solid line). 

of only 0.015 for the dust component and of 0.08 for the 
synchrotron component, exploiting polarization data only 
(Planck plus the WMAP 23 GHz channels). 

Perhaps the most important new result of this paper is 
the elaboration, for the first time, of successful methods for 
estimating component separation errors. We presented two 
alternative error estimation methods for the CCA, relying on 
completely different assumptions (one is based on studying 
the marginal probability for each parameter estimated, the 
other on the redundancy of results for different patches of 
the sky) so that they can also bo used to cross-check the 
results. A first application of these methods allowed us to 
identify the most uncertain values of the parameters and to 
compute the appropriate weights to be used to combine the 
estimates in different regions of the sky to build a smooth 
map of foreground spectral indices. 

The components were reconstructed with a Generalized 
Least Square (GLS) solution in pixel space, which allowed us 
to fully exploit the spatially-varying information obtained. 
Even if the choice of the channels used in the reconstruction 
is optimized for the CMB, wc were able to reconstruct a 
very accurate dust polarization map (correlation coefficient 
~ 1 with the input maps). Errors in the estimation of the 
foreground spectral indices were successfully propagated to 
the foreground maps. They turned out to be well below those 
due to instrumental noise, except for the dust component 
close to the Galactic plane. 

As for the CMB, we used ROM Aster to estimate the 
polarization E-mode power spectrum from the reconstructed 
CMB map and found negligible effects by the residual fore- 
ground contamination even masking only 10% of the sky. 
Again the component separation errors, propagated to the 
power spectrum, were found to be subdominant with respect 
to noise. We also showed that our component separation 
method can be useful to tackle B-mode detection, once the 
experimental noise level allows it. 
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